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Abstract 

The exact computation of the nearest-neighbor spacing distribu- 
tion P(s) is performed for a rectangular billiard with point-like scat- 
terer inside for periodic and Dirichlet boundary conditions and it is 
demonstrated that when s — > oo this function decreases exponentially. 
Together with the results of Ref. [13] it proves that spectral statistics 



of such systems is of intermediate type characterized by level repul- 
sion at small distances and exponential fall-off of the nearest-neighbor 
distribution at large distances. The calculation of the n-th nearest- 
neighbor spacing distribution P n (s) and its asymptotics is performed 
as well for any boundary conditions. 



Unite Mixte de Recherche de l'Universite Paris XI et du CNRS (UMR 8626) 



1 



1 Introduction 



The statistical analysis of quantum energy levels for a given system in the 
semi-classical limit is a well-studied feature in the theory of spectral statistics 
[0]"@- The main conjectures in this field are the follow: 

(i) The Berry- Tabor conjecture Q: generic integrable systems obey Pois- 
son statistics, which implies that their energy levels behave as independent 
random variables. 

(ii) The Bohigas-Giannoni-Schmit conjecture ||: generic chaotic systems 
follow the Wigner-Dyson distributions of Random Matrix Theory (see ||). 

There is enormous amount of numerical evidences that many physical 
systems do agree with these two main level statistics. Partial analytical 
results support these conjectures for integrable rectangular billiards || and 



quantum chaotic systems [M- 10 



However, there exist systems which are neither integrable nor chaotic and 
their spectral statistics do not follow any of the above leading models. In 
many cases their statistics have features intermediate between the Poisson 
statistics and that of random matrix ensembles and for this reason they are 
called "intermediate statistics" |TT| , fT2|l , [|13|]. For the first time this type of 
behavior had clearly been observed numerically for the 3-dimensional Ander- 
son model at the metal-isolator transition point [[ll|] and later it was argued 
bf] that spectral statistics of diffractive and pseudo-integrable systems is 



also of intermediate type. 

To study precisely the statistical behavior of the energy levels of quan- 
tum systems one usually introduces different functions that characterize the 
statistics . The most important quantity for our purpose is the distribution 
of nearest-neighbor spacings, -P(s), which is the probability that two levels 
are separated by a distance s with no levels inside this interval. 

For the Poisson statistics the nearest-neighbor distribution takes the fol- 
lowing particularly simple form 

P(s) = exp(-s) (1) 

and it is characterized by (i) the absence of level repulsion (-P(O) ^ 0) and 
(ii) exponential decay for large distance. 

For standard random matrix ensembles the nearest-neighbor spacing dis- 
tributions are given by complicated expressions Q but their main features 
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can be seen from the Wigner surmise 



P( 



aps^ exp(— cps' 



(2) 



where j3 — 1,2 and 4 corresponds, respectively, to orthogonal, unitary or 
symplectic ensembles, and ap and cp are constants determined by the nor- 
malization conditions. Its main properties are (i) level repulsion, P(0) = 0, 
and (ii) very quick decrease at large values of s, P(s) oc exp(— cs 2 ) when 
s — > oo. 

We call spectral statistics of intermediate type if they have the follow- 
ing hybrid properties (cf. [11], Q]): (i) they exhibit the level repulsion, 
P(0) = 0, as for standard random matrix ensembles, and (ii) they have ex- 
ponential decay at large s, P{s) oc exp(— cs) when s — > oo, similarly to the 
Poisson statistics. Little is known analytically for systems with intermediate 
statistics though it is possible to write down models which will have this type 
of statistics |i~5| , Hl2| . 



The rectangular billiard with a point-like scatterer inside belongs to the 
class of diffractive systems and is one of the models which is supposed to 
have intermediate statistics ill], [I3|. Without the scatterer this model is an 
integrable system and when the ratio a 2 /b 2 of the sides of the rectangle is 
a 'good' irrational number its quantum energy levels {e n } obey the Poisson 
statistics M, 101. The addition of a S— function scatterer 



V = X S(x — xq) 



(3) 



inside the rectangle corresponds to a rank-one perturbation, and the new 
quantum energy levels E of the perturbed rectangular billiard will obey the 
following quantization condition |TB[, [TTJ 



(4) 



where ipn^ and e n are the eigenfunctions and eigenvalues of the unperturbed 
rectangular billiard. Similar equations appear in different models. In partic- 
ular, the quantization condition for the star graphs can be transformed to 
this form [18 
as well. 



All our results are applicable without changes in such cases 
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Rank-one perturbations have been studied in the context of ballistic mo- 



tion of particles in regular [17] or chaotic cavities [19], and in the context of 



random matrix theory PD| . When a 5— function potential (^) is added to a 



chaotic system with random matrix statistics it has been proved pi that the 



new eigenvalues in Eq. (J3]) are also distributed according to the same statis- 
tics. In the chaotic case the universal part of the spectral statistics is not 
changed by finite-rank perturbation. On the contrary, when the unperturbed 
system is integrable, the perturbation (|j) changes dramatically its spectral 
statistics |L7], |jl3 . 



In Ref. [[13[] the two-point correlation function of a rectangular billiard 
with a small-size scatterer inside (described by the potential (Q)) has been 
computed analytically. One of the conclusions of this paper was that spectral 
statistics of such singular billiards do have level repulsion. For billiards with 
periodic boundary conditions the two-point correlation function and, conse- 
quently, the nearest-neighbor distribution vanish linearly at small distances 
with the slope independent on the coupling constant 

P{s) -> —s. (5) 

For billiards with Dirichlet boundary conditions the small-s behavior of the 
two-point correlation function and the nearest-neighbor distribution is dif- 
ferent: one has (see |T3| ) 

p( s )n°-L s log 4 S . (6) 

The main purpose of this paper is to compute analytically the nearest- 
neighbor spacing distributions for this model and demonstrate that any 
boundary conditions they decrease exponentially at large separation. To- 
gether with the results of Ref. |1| about the existence of level repulsion it 
will furnish the proof that spectral statistics of these systems are of interme- 
diate type. 

The plan of the paper is following. In Section |2] we generalize the for- 
malism used in [|1^] to describe the nearest-neighbor spacing distribution 
for a billiard with a point-like scatterer for periodic and Dirichlet bound- 
ary conditions. Though the resulting formulas are explicit and exact they 
are quite cumbersome and in Section ||] we study the asymptotic behavior of 
P(s) for large s. It is demonstrated that in all cases the nearest-neighbor 
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distribution has exponential tail at large distances thus proving the inter- 
mediate character of spectral statistics of singular billiards. In Section ^ 
the n-th nearest-neighbor spacing distributions for these billiards with peri- 
odic and Dirichlet boundary conditions are computed analytically and their 
large distance asymptotics are found as well. In Appendix we present certain 
technical details of the computation of necessary integrals. 

2 The general formalism 
2.1 Preliminary computations 

In this section, our aim is to find analytical expressions for the nearest- 
neighbor spacing distribution of the solutions, E, of the following equation 

JV 

AV-^=1 (7) 

where ej j = 1, . . . , N are independent random variables with a uniform 
distribution dfi(e): 

da(e) = l^ de *-W<e<W 
| otherwise 

and Tj are positive constants with mean value 1 

1 N 

n=l 

This normalization condition permits to introduce conveniently the coupling 
constant A. 

In general this equation describes zeros of a meromorphic function whose 
poles are assumed to be independent random variables and it can correspond 
to different physical problems (see e.g. f22fl). In this paper we prefer to 
consider it as the quantization condition of rectangular (or more general 
integrable) billiards with a small-size impurity inside ||16|| , [|17]]. To ensure 
that energy levels of unperturbed billiards behave as independent random 
variables it is necessary to assume that the ratio of squares of the sides of the 
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rectangle, a and b, is an irrational number badly approximated by a rational 
(that is a diophantine number) || with the following property 

, a 2 m, C , , 

> 3 ( 10 ) 



b 2 



n n n 



for all integers m, n and some k > 2. 

The residues, r n , depend on boundary conditions. For quantum problems 
with periodic boundary conditions r n = 1. For Dirichlet conditions 

a ■ 2 n ■ 2 71 rtt\ 

r mn = 4sm -mx sm -my (11) 
a b 

where xo, yo are coordinates of the singular scatterer. When the ratios xo/a 
and yo/b are non-commensurable irrational numbers and m, n — > oo r mn can 
be considered as independent random variables 

r mn = 4 sin 2 0i sin 2 02 (12) 

with angles fa uniformly distributed between and tt/2. 
When both ratios x^/a and yo/b are rational numbers 

xo _ Pi _ P2 (13 , 

a q\ b q 2 

with co-prime integers (pj, g,) the residues (|ll|) only depend on m mod q± 
and n mod q 2 and there is only finite number of residues determined by qi 
angles fa in fljjD (see |13] for more detail) 



ki 

fa = tt — , with ki — 0, 1, . . . ,qx — 1; 

^2 

02 = 7T — , With A; 2 = 0, 1, . . . ,g 2 — 1- (14) 
92 

All our formulas below remain valid for general r n . 

Obviously there are iV solutions Ej of Eq. (|7|) since each interval ]e$, e^i[ 
contains one and only one of these solutions. We are interested in the nearest- 
neighbor distribution, P(s), that is the probability that two energy levels 
Ei and Ej are neighbors separated by a distance s. In our case it is the 
probability that 2 solutions Ei and Ej of (0) are separated by one and only 
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one unperturbed level e&, and that \Ei — Ej\ = s. Let us compute at first the 
probability P(Ex,E 2 ) that two given energy levels E\ and E2 be neighbors. 
Assuming for instance that E 2 < Ex, P(Ei,E 2 ) is the probability that that 
one solution of (|7|) equals Ex, another one equals E 2 , and that there exists i, 
1 < i < N, such that 

eie]E 2 ,Ex[ , v 

\/j^i,e j ^)E 2 ,E 1 [. [i0> 

As it is supposed that are independent random variables with a uniform 
distribution, 

rW N 1 N 

P(Ex, E 2 )= J] ^p(E 1 )p(E 2 ) £ X(ei) R (1 - x(e 3 )) , (16) 

fc=l i=l j^i 

where xi e ) is the characteristic function of the interval }E 2 , Ex[ equal to 1 if 
e belongs to ]E 2 , Ex[, and to otherwise, and p is the density of the solutions 
En 

N 

p(E) = J2$(E-E l ). (17) 

i=l 

It is convenient to rewrite these formulas in a more symmetric way: 

N 

P(Ex,E 2 ) = J2 \{d^ k {ek)p{E l )p{E 2 ), (18) 

where variables a^, k — 1, . . . , N take 2 values: 0, 1; and we introduce two 
different measures 

d^o(e)0(e) = — y ^x(e)0(e)rfe = — <f>(e)de (19) 



and 



W 1 /-W /-Bi 



d//i(e)0(e) = — / (1 - x(e))0(e)cfe = — ( / - / )0(e)de. (20) 
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The summation in flIB|) is performed over all sequences at which contain 1 
zero and N — 1 ones. 

Because {Ei} are solutions of Equation (|7|) the density of states ( |T7| ) can 
be rewritten under the form (cf. fl3| ) 



N ,\N 



Representing the 5— function by a Fourier integral one gets 



--- - - \ \ i=i 



271 

and finally the probability (|T6[) can be put under the form 



x 

k 



P( El ,E 2 ) = r d -^J2 [ fld.Me-^^ (23) 

J —CO ^ r , J i i 

Wj fc— 1 

iV AT , 



, AEi - e kl ) 2 (E 2 - e fc2 ) 2 . 11 \ E x -e jx E 2 - e h 

Let us introduce the following functions: 

/ -(ai,a 2 ) = dfi a (e) exp ( z— — ^ h £- 



* iff (ax,a a ) = /^(^^exp^-^ + i-^j, (24) 

^ (ttl ' aa) = / ^ (6) (^ - e)'(E 2 - e )=» 6XP ('E^ + % ET^) ■ 

The nearest-neighbor distribution can be expressed through these functions 
in the following way 

p( El , e 2 ) = y: r ^ ( e ^ M n to) ( 2s ) 

+ E r ^i-> r ^2 CTfc K fc ) n A to) ) e" i(ai+a2)/A . 
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Here and below when it will not lead to a confusion we use the notation f(a) 
for a function of two variables f(ax,a 2 ) and /(or) instead of f(a±r,a 2 r). 
This formula is valid for all sequences of a^. The functions ( ^ip with index 
a = correspond to unperturbed level between E 2 and the functions with 
index o = 1 correspond to unperturbed level outside ]E 2 ,Ei[. Therefore to 
describe the nearest-neighbor distribution the summation should be done 
over iV possible sequences containing only 1 zero. 

The functions and g a are related to f a by the relations 

d 2 <9 4 
^ ja (a 1 ,a 2 ) = -^2/<r(ai,a 2 ), g<r{ai,a 2 ) = g 2g 2 ^( tt i' Q2 )- ( 26 ) 
i 12 

Therefore in order to compute P{E%, E 2 ) it is necessary to find only f a . Let 
us introduce the functions 

I (T (a 1 ,a 2 ) = 2W [ d Ay ( e )[l- exp (i_^_ + i_^_)], (27) 
which are related to our basic functions f a as follows: 

= 1 -2^"2W /l(ai ' a2) ' 
/o(ai,a 2 ) = — / (tti,a 2 ), (28) 

where u; = .Ei — .E^ is the difference of energies (we recall that we have 
assumed E 2 < Ei). 

The integral defining I\(a) can be split into two parts 

f E2 f W "i « 2 

7i(ai,a 2 ) = (/ +/ )[l-exp(z— h i— )}de 

J-w J Ei E x -e E 2 - e 

= Ji(a ll a 3 )+j(a 1 ,a 2 ) 1 (29) 



where 



and 



f E * f°° ai a? 

Ji(a 1 ,a 2 ) = ( +/ )[l-exp(i— — + i-^ —)]de (30) 
J-oo Je, E x - e E 2 -e 



/—W poo 
+ )[l-exp(i— ^ + i—^-)]de. (31) 
Jw Ex-e E 2 -e 
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For convenience we define the function Jo (a) = Io(a) so that from (E 



f El «i a 2 

J ( ai ,a 2 )= [l-exp(z— )}de. (32) 

Je 2 Ei\ — e E 2 — e 

The integral (|3l|) defining j(a) has no singularity inside the integration region 
and as it was demonstrated in |I3| it is sufficient to take into account only 
terms linear in a and to ignore the difference between E\ and E 2 (i.e. set 
Ei w E% w E). In this approximation 

W - E 

j(a u a 2 ) =i( ai + a 2 )\og w - — . (33) 

On the contrary the functions J a (a) are quite cumbersome. One can easily 
check that they depend only on the difference of energies, u = E% — E 2 , and 
that 

rv 

J a (a) = uJ a (-) (34) 



where the functions J a (a) are defined by Eqs. (p0|) and (|32|) with E\ = 1 and 
£ 2 = 0. 

In Appendix it is demonstrated that these functions obey the differential 
equation 

(dt - d 2 )J a ( ai , a 2 ) = e i(Ql - a2) 0<x(«i, « 2 ) (35) 

where <9j denotes the derivative with respect to and functions CT (a) at real 
a are given by ( |103j ) and ( |106[ ). From this equation it follows (see Appendix 
for details) that the function J\{a) is an analytical function of 2 complex 
variables ati, a 2 with the cuts as in Fig. |l|a and [l|b given by the following 
expression 



dt 



a 2 ) = j ^ (1 - e l(Ql+a2)t ) + / e i(2i - ai " a2 Vi(t, «r + « 2 - t)tft, 

(36) 

where J- denotes the principal part of the integral. 



10 



a 



«2 



(a) 



(b) 



a, 



a 2 



(c) 



(d) 



Figure 1: The cuts in complex planes of a.\ and a 2 for functions Ji and 4>x- 
(a), (b); J and (p : (c), (d). 



The function Jo(a) is an analytical function in a region indicated in Fig. |TJc 
and [l|d with integral representation 



JoOi,a 2 ) 



)+/ * 





< 2 *- ai - a2 Vo(i, ai + a 2 - t)dt. 

(37) 



Exactly as it was done in [O] one can prove that functions ^i a , g a defined 
in Eqs. ([24]) can be expressed through the functions <p a as follows: 



g a {a 1 ,a 2 ) 
*i <T (a 1 ,a 2 ) 
*2<j(a; 1 ,a2) 



a; 



1 ,< 



2iy 



.«! a 2% 

• u7' 77' 



(38) 



In, /"l «2- 
^20a — , — , 
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2.2 Nearest-neighbor spacing distribution 

Using Eqs. ( p8[ ) one can integrate the first term in ( p5|) by parts and because 
(see (H) and (H) 

(ft - d 2 )Ua u a 2 ) = -J-e^- a ^U-, -) (39) 

Z Co 1 Cl> 



one obtains 



r/nif/a2 E^(^)[n^(^)] e " i(Qi+ " !iA 



47T 



2 



°° """ 3 k^j 

oo 



x [nAK')] e ~ i(Qi+a2)/A - ( 4 °) 



According to Eqs. (p58|) the second term term in can also be expressed 
through the same functions </> CT (a) and after the scaling of variables a — > au 
(cf. (0)) the nearest-neighbor distribution (|25|) takes the form 

x [n/,h)]e-" (ttl+ ° 2)/A , (41) 



where /^(a) = f ai (au>). Using Eqs. (pq), (p9|), and (|31|) one gets 

u; W - E 

/i(ai,a 2 ) = 1 - —(1 + Ji(ai, a 2 ) + i(ai + a 2 ) log — - ^ ), 

/o(ai,a 2 ) = 2^(1- Jo(ai,a 2 ))- (42) 

Expression ([H|) is valid for any sequence {<Tk} G {0, 1}^. To get the nearest- 
neighbor distribution one has to sum over N sequences containing only one 
zero. Taking into account that in the limit N — > oo the restriction j ^ k is 
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unessential we obtain, keeping only the dominant term, 
N 2 f°° da.\da.2 ~ 



X 

1 



[]/i(an)]e— (Q1+Q2)/A , (43) 



where the operator Vf )B (a) is defined for arbitrary functions f(a) and g(a) 
by the following expression 

V f>g (a) = < rf(ar)e i{ai - a2)r >< rg(ar)e i{ai - a2)r > 

- < (— /(or))e^^> >< (j— g{ar))e^- a2 > >, (44) 

OOtl 00(2 

and < f(r) > means the mean value over all values of r 

1 - 

</(*•)>= ^E/W' ( 45 ) 

n=l 

Measuring the energy difference u in the units of mean level spacing 

Nuj 

s = (46) 

2W V ; 

the product Yli fi( ar i) can also be simplified in the limit of large N (see 

T^r ~ iVcu - W - E 

11 /i(arO « ^(-^yt 1 + J i(oi» «2) + i(ai + a 2 ) log — — ]). (47) 

Introducing the renormalized coupling constant A' 

— = hog (48 

A' iVA 1¥ + £ 1 ; 

we obtain the final formula for the nearest-neighbor distribution P(s) = 
{2W/N) 2 P{uj): 

+V*,,, 1 (a) + ^^ («)] e - s « Jl ^ >+4 ( Ql+Q2 )/ A '). (49) 
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2.3 Analytical continuation 

Usually, if one wants to compute an integral 

/oo 
da 1 da 2 f(a 1 ,a 2 )e~ sJ{aua2 \ (50) 
-oo 

where J (a) and /(a) are analytical functions in a certain region, the first step 
is to move the integration contour as far as possible to decrease the integrand. 
In general, during that deformation one can either meet a saddle point or 
a singularity which signifies that further deformation of the contour either 
will increase the integrand or is not possible. If no such obstacle appears the 
integral is zero. 

In the case of the nearest-neighbor distribution (^9|) the saddle-point 
equation reads (taking here r = 1) 

~ i d ~ i 

Ji(«i, a 2 ) + — = 0, - — Ji(«i, a 2 ) + — = 0. (51) 



da i A' da 2 A 

In particular, these equations imply that at any saddle-point 

, d d 



1 dazi da 2 



)J 1 (a u a 2 ) = 0. (52) 



From ( p. 00] ) and ( |102j ) it follows that this difference is proportional to Kq(x) 



with x = 2y/aia 2 . But K (x) has no zero on the complex plane (see ^3] 
p. 62). Therefore our integral (|49|) has no saddle-points and one can move 
the contour of integration freely. If the pre-factor in (|50) has no singularities 
in the region where J\ is analytical the contribution vanishes. Note that 
it is exactly what should be expected for physical reasons: replacing the 
pre-factor in ([K]) by 1 we have to obtain the probability that there are two 
perturbed levels in E 2 and E\ but no unperturbed energy levels between, 
which according to Eq. (^) is impossible. These considerations prove that 
the integral ([49|) with a pre-factor equal to 1 or, more generally, with any pre- 
factor analytical in the same domain as J\{a) (and not too quickly increasing 
on infinity) must vanish. 

These arguments allow to simplify considerably the expression ( [4"5D for 
the nearest-neighbor distribution. The pre-factor in Eq. ([|£]) is 

f(a) = s(l- < J (ar) >)V r 01i01 (a) + V Ml {a) + V^ (ct). (53) 
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First, as only the functions with index have analytical properties different 
from that of J\{at) (see Fig. p, one can drop the first term and keep only 

f(a) = -s < J (ar) > V^,^(a) + ^^(a) + V^ >lf>0 (a). (54) 

Secondly, according to Eqs. ( |107|) and ( |117| ) 

4> (a) = -0i (a) + 7r(sgn(ai) - sgn(a 2 ))io(2v / -"i«2), (55) 

and 

Jo(a) = —J\{a) + sgn(«i)-R(a) + sgn(a 2 )i?^(a), (56) 

where the function R(a) is defined in ( |118|) and R^a) = R*{a 2) a\). 

When these expressions are substituted in (|54]) the terms with index 1 
can be dropped out because they have the same analytical properties as 
Ji(a) and, as it has been discussed above, their integrals vanish. Finally the 
pre-factor takes the form 

f(a) = -s [sgn(ai) < R(ar) > +sgn(a 2 ) < R\ar) >] (a) 

+ 7r(sgn(cn) - sgn(a 2 ))(V r j 0i ^ 1 (a) + V^ 1>Jo (a)) (57) 

+ 2m^-(< e** >) 2 - 2m^-(< e^ r >) 2 , 
a 2 + it oil — it 

where Jo is the Bessel function Jq(2^ —aia 2 ) . The last term in this equation 
appears from the differentiation of (sgn(ai) — sgn(a 2 )) in ([44]) . The resulting 
5-function allows us to take the remaining terms at small a. We also write 
«i,2 T * e where e — > + to remind the region where the functions are defined. 

When this expression is substituted in Eq. P9| ) one gets the final formula 
for the nearest-neighbor distribution which can be conveniently written as a 
sum of 3 terms 

P(s) = ^- 2 (A(s) + B(s) + C(s)), (58) 

where 

A(s) = -s J °° dan da 2 < R(ar) > V Ml {a)e-< <J ^ ra > +i ^ +a ^l x '^ 
-sJ™da 2 JZ da 1 < R){ar) > V Ml {a)e-« J ^ ra > +l ^ +a ^l^ + c.c. (,59) 
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B(s) = vr da, da 2 (V^ Jo (a) + Vj ^(a))e^< J ^> + ^ + ^ 
-7T da 2 &i(^,j„(«) + K7 Oi01 (a))e-« J ^ ra )>^^ + ^)/ A ') + c.c.(pO) 
and 



OO J 

"" 2 fV p -ia 2 r ^\2_-s(<Ji(0,ro 2 )>+ta 2 /A') 



C(s) = 2tu / — (< e- M2r >) 2 e 

-oo «2 + «e 

r"0O 



_ 27rz /°° _^^ (< e ia ir > N2 e - S «J 1 (m lj 0)>+ia 1 /V)_ (61) 
J-00 "1 - « 

These expressions look quite complicated but in the next Section we show 
that their asymptotics when s —>■ oo can easily be computed. 

3 Asymptotic behavior 

The formulas (|58|)-(|6T|) have been written in such a way that when the inte- 
gration is performed from — oo to +00 the contour of integration as a whole 
can be shifted into the complex plane. The direction of such a deformation is 
different for the integration over a, and that over a 2 . It can conveniently be 
fixed by the following change of variables: a, = —iv or a 2 = iv. In the new 
variable the allowed deformation of the contour is in both cases Rev > 0. 

Let us consider first the simple integral ([H]). From Eqs. ( |113| ) and ( |119| ) 
of Appendix it follows that Ji(0, iv) = Ji(—iv,0) = I{v) where 

f°° dt 

I(v) = i (l- e - vt ) = e v -l-vEi(v), (62) 



and 



1 t 2 



Ei(v) = -fje-t (63) 



is the standard exponential integral (see e.g. p3| p. 143). 

Consequently, after the above change of variables the integral (^T|) takes 
the form 



C(s) = 


-2m 


p+lOO 


dv 


;< e vr >)V 


-s(<I(vr)>-v/\') 




' —ioo 


v 






-2ni 


p+ioo 
' —ioo 


dv 
v 


;< e vr >)V 


-s(<I(vr)>+v/X') 



(64) 
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Now one can move the contour to the right till it goes through the saddle 
point. For the first integral the position of the saddle point, v+, is defined 
by the equation 

Tv < /(t '+ r) >= \> (65) 

and for the second integral the saddle point v _ is determined from the similar 
equation but with changed sign of the coupling constant A' 

Tv < ,{vr) >= "F (66) 

As I'(v) = — Ei(v) the saddle points v± are roots of the equation 

< rEX(v±r) >= ^— (67) 
A 

and it is possible to prove that for any real A' there is one and only one 
solution of this equation. 

Expanding the exponent in (|6^ ) in vicinity of the saddle points and taking 
into account that I"{v) = —e v /v and < I(rv±) > =Ft>±/A' =< e rv± > —1 one 
gets that in the limit s — > oo the function (|61|) is the sum of contributions 
from 2 saddle points 

C(a) = E u "J e-«<^>-» (68) 

f—/ J< re rVi > svi 

The saddle point with the smallest value of < e rv >, which we denote by v sp , 
dominates and it should formally be the only one to be taken into account. 
For finite A' it corresponds to the solution of Eq. ( |57|) with negative right- 
hand side 

< rEi{rv sp ) >= -— (69) 

Note the appearance of the absolute value of the renormalized coupling con- 
stant. When A' — > oo both saddle points v± will give comparable contribu- 
tions and both should be included. 



The asymptotics of the other terms ([59]) and ( p0|) can be computed by 
similar considerations. These functions are defined as double integrals, the 
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first one is taken from to 00 and the second from —00 to 00. To compute 
their asymptotic behavior for large s the latter integral should be deformed 
into the complex plane as it was done above and in the former integral one 
has to take into account only the lowest order terms according to expansion 

Let us consider the first term in Eq. ( |59|) . We need to know the limiting 
behavior of the integrand when a\ — > and ct2 + ol\ = iv (we prefer to use 
this deformation instead of the usual one, «2 = iv, to simplify the formulas 
below). From Eqs. ( |113| ), ( |122| ), and (|126|) it follows that for real a 

Ji(a x ,iv - «i) ai ^° I(v) + e v [^- - i{2^ + \ogv - 1 + logai)]«i, (70) 

with I(v) defined by (p2[). From ( |118| ) in this limit R(a) — irotie v . The 
dominant contribution in V^ 1)( ^, 1 (a) comes from the second term in ([0]) and 
one gets V^^^a) ~ — i < e vr > 2 /(a%v). Combining all terms together and 
changing the variable a\ — > a/(s < re rv >) we find that at large s 

poo r+ioo 7 

A(s) = 2m / da / — (< e™ >)2 e -«/(«-)WA') 

JO J-ioo v 

xe-™ /2 sin[a(loga-logs + #(i;))] + (A' — > -A'), (71) 



g(v) = logv + 2 7 - 1 + 2 B - log < re™ > . (72) 

< re'" > 



where 



The integral over v is an analog of Eq. ( |64l) and can be computed exactly as 
above: 

A(s) = -(2vr) 3 / 2 ^^^L / (log S -,K))e-^>- 1 ), (73) 
where the function f(y) is given by the integral 

POO 

f(y) = / e^ 2 sin[«(loga - y)]da. (74) 
Jo 



One can check that the contribution ( 160]) when s — > 00 is smaller by a factor 
1/s with respect to fl6"8|) and (fT3|) and can be neglected. 



18 



Finally, we obtain that the nearest-neighbor distribution, P(s), in the 
limit of large s has the following asymptotics 

( < p rv sp > \2 p -s<e rv "v > 

P{s) = j== = } = : t=— [1 - /(logs - 9(v sp ))]. (75) 

The saddle point value, v sp , depends on the renormalized coupling constant 
by Eq. @ 

- < rEi(rv sp ) >= — (76) 

When A is very large the contribution of the second saddle point with reversed 
sign of right-hand side of this equation should be added. 

In Fig. ^]the plot of the function f(y) defined in (|7|) is presented. When 
y — > oo this function goes to zero as —1/y. Therefore the true asymptotics 
of P(s) is given by the first term in ( [75] ) 

P(s) = / [<C >) =- (77) 

But because in Eq. (|75|) the argument of the function f(y) up to the constant 
( |72|) is logs this decrease is quite slow and at numerically accessible values 
of s of the order of 10 (i.e. y of the order of 3-4) as it is evident from Fig. |2] 
this function gives a noticeable contribution. 

For the billiard with periodic boundary conditions the residues in (|7|) all 
equal 1, and all mean values are reduced to the corresponding function, i.e. 
for any function / the quantity < f(rx) > becomes f(x). 

In this case the nearest-neighbor distribution has the following asymp- 
totics 



e 3v ap e -se^P 



P{s) = ^2^~^/T [l - f{logs ~ 9{Vsp))] - (78) 

The value of v 8p is determined by the equation 

HO = -i^i . (79) 
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Figure 2: Plot of the function f(y) denned in (T73). 



and 



9\v, 



spy 



log V, 



sp 



' sp 



2 7 - 1. 



(80) 



In Fig. [| we present the comparison between numerical computations 
and the theoretical prediction ( |78|) for 2 values of A. The logarithm of the 
nearest-neighbor distribution is plotted as a function of s. The upper curve 
(squares) corresponds to 10 6 levels with A' = 1 and the lower one (circles) to 
5.10 5 levels with A' = 100. For clarity the curve corresponding to A' = 1 has 
been shifted on the right by s — > s + 2. The solid lines represent theoretical 
predictions fl78|) for these values of the coupling constant. For A' = 100 
two saddle points with different sign of the right-hand side of Eq. ([T^) have 
been taken into account which roughly doubles the result ( |78| ) . These results 
very well confirm theoretical asymptotics of the nearest-neighbor distribution 

In the case of a rectangular billiard of size a x b with Dirichlet boundary 
conditions, with a point-like scatterer such that the ratios of its positions 
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2 4 6 8 10 12 

s 



Figure 3: The nearest-neighbor distribution in the periodic case. Squares 
and circles correspond respectively to A' = 1 and A' = 100. Solid lines: 
theoretical predictions (f78| ). For clarity the upper curves are shifted to the 
right by 2 units. 



(xq, Do) to the corresponding sides are non-commensurable irrational numbers 
the residues, r n , can be considered as random variables of the form given by 
Eq. ([12] ) and the mean value of a given function /, < f{r) >, should be 
computed as follows 

< f(r) >=— 2 \ l d<\> x \ l d<f) 2 /(4sin 2 lS in 2 2 ). (81) 
71 Jo Jo 

In the case where x^/a and yo/b are rational numbers (13), the mean value 
< f(r) > takes the form (see fTB]j ) 



91— 1 92-1 



< f(r) >= — J2H /( 4sin2 — sin2 — )■ ( 82 ) 



Q1Q2 — — Qi Q2 
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The cases ki = and k% = correspond to unperturbed subsequence of levels 
(for Dirichlet boundary conditions) and one has a freedom to include them 
to the level density ([TJ]) or consider them separatively. In the later case the 
terms with ki = should be omitted and l/^iffe) in front of the sum should 
be substituted by l/((gi — l)(g 2 — l)). The generalizations of < f(r) > for the 
case when only one ratio Xo/a or y Q /b is an irrational number and the other 
one is a rational number or both ratios are irrational but commensurable 
numbers are straightforward. 

With such definition of the mean value the asymptotics of the nearest- 
neighbor distribution is given by Eq. (|75|). 

As in the periodic case the quantity — log P{s) is expected to be nearly 
linear with the slope < e™ sp > where v sp is the solution of Eq. (|TB|) . 

Fig. H shows the comparison between numerically computed nearest- 
neighbor distribution for a billiard with Dirichlet boundary conditions and 
the expected asymptotic behavior (|75| ) (solid line) for 2 values of the renor- 
malized coupling constant, A' = 1 and A' = 10 (xo/a, yo/b are incommen- 
surable irrational numbers). To better understand the asymptotics of the 
nearest-neighbor distribution we present in Fig. [5] the functions which de- 
termine the exponential decrease of P{s). Different curves in this figure 
correspond to the following functions: — Ei(x) (thick solid line), expx (thin 
solid line), — < rEi(rx) > (thick dashed line), and < exprx > (thin dashed 
line). 

For periodic boundary conditions the intersection of the horizontal line 
having the ordinate 1/A' with the graph of — Ei(x) gives the value of v sp . The 
point of intersection of the vertical line going through v sp with the graph of 
e x determines the exponent in Eq. (|78|). For Dirichlet boundary conditions 
one should use the same procedure but with the functions — < rEi(rx) > 
and < e rx >. 



4 The n-th nearest-neighbor spacing distri- 
bution 

In the previous Sections we compute the nearest-neighbor distribution, P(s), 
for a rectangular billiard with a small-size scattering center inside, i.e. the 
probability that 2 levels are separated by a distance s with no levels in 
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11 13 



Figure 4: The nearest-neighbor distribution for a billiard with Dirichlet 
boundary conditions. Squares and circles correspond respectively to A' = 1 
and A' = 10. Solid lines: theoretical asymptotics flT8|). Upper curves shifted 
to the right by 2 units. 

between. In this Section we generalize the formalism to compute the n-th 
nearest-neighbor distribution, P n (s), which is defined as the probability that 
2 levels at distance s are separated by n levels, for n > 1. 

Our starting point is the expression (^TJ) which is valid for any sequence 
<Jk of and 1. To obtain the n-th nearest-neighbor distribution one should 
sum over all sequences of length N with exactly n + 1 zeros. Performing the 



same steps as in Section 2.2 and taking into account that when N — > oo we 
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Figure 5: Functions — Ei(x) (thick solid line), — < rEi(rrr) > (thick dashed 
line), e x (thin solid line), and < e rx > (thin dashed line). 

have C%j — > N n /n\ one obtains 

PM = e L < Mar) >] v ** (a) 

s n 

+—(1- < J (ar) >) n (V Ml (a) + V Mo {a)) 

+7^tt(1- < JoM >)»- 1 ^^(a)]e-« Jl ^» +i ( ai+aa )/ v ), (83) 

with all functions defined as above. 

As a consistency check one can verify that the sum over all n coincides 
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with the exact expression of the 2-point correlation function derived in 



R 2 (s) = Y,Pn(s). (84) 

n=0 

Similarly to Section |2.3| the different terms in Eq. ( |33D can be classified 
according to their analytical properties. Three groups of terms appear. First, 
terms which have the same analytical properties as the function J\{ot) in the 
exponent. Exactly as it was done in Section |2.3| one can argue that their 
contribution is zero because one can freely move the integration contour to 
infinity. The second group consists of terms which are singular on one variable 
but have 'good' analytical properties on the other variable of integration. 
The large s asymptotics of such terms can be calculated as in Section |3] 
by computing the leading terms over the former variable and shifting the 
contour of integration over the latter variable into the complex plane till it 
reaches the saddle point. The asymptotics of these terms will be proportional 
to exp[— se Vsp ] as for the nearest-neighbor distribution. Finally, the third 
group (which does not exist for P(s)) includes terms which are singular (i.e. 
have analytical properties different from that of Ji(oe)) for both variables of 
integration. In this case no deformation of the integration contour is possible, 
the region of small ot\ and a 2 will be important and the asymptotic result 
will be proportional to exp(— s). As v sp > it is the last group of terms 
which will dominate the asymptotics of P n (s) when s — > oo. 



From formula (122) of Appendix it follows that, when both variables a>i 



and a 2 are small, the functions V^ it ^ will be equivalent to their singular parts 

1 2ni 

VW,*i( a ) ' V Mo ( a ) S(a 2 ), 

ol\Ol 2 o.\ 

2iri 

V^,0 o (a) ~ 5(ai), V Mo (a) ~ (2n) 2 S(a 1 )S(a 2 ). (85) 

a 2 

The singular terms in (1— < Jo(ra) >) n are given by the small a behavior 
Eq. (121). Keeping only the terms 1 
^< J (ra) > + < Ji(ra) >)] n we get 

1- < j (rat) >) n = 1 - n7r(|ai| + |ot 2 |) + n(n - l)7r 2 |ai||a 2 |. (86) 



in Eq. (121). Keeping only the terms with a singularity in [1+ < Ji(ra) > 
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(We recall that the terms coming from J\ vanish.) The first term in this 
expression is the dominant regular contribution, the second one is the domi- 
nant contribution singular in one variable, and the third term is the dominant 
contribution singular in both variables. 

Combining all the above expressions together and using (|120 ) we obtain 

s n ~ x _ s f°° da 1 da 2 r 2 2 / \ / \ 
P n (s) = J ^ -^[vr s sgn(a 1 )sgn(a 2 ) 

+27T is sgn(ai )S(a 2 ) - 2it 2 is sgn(a 2 )5(a 1 ) + {27r) 2 5{a 1 )5{a 2 )} (87) 



7T 

x exp -s[-(|«i| + \a 2 \) - iax(log \ai\ + + ia 2 (\og \a 2 \ + g + )\, 



where 



g± =< r logr > +7 - 1 ± — . (88) 

A 

After simple calculations we obtain the large s asymptotics of the n-th 
nearest-neighbor distribution 

„n-l 

Pn(s) = 7 rrre-'fl - /(bg S ~ 9+ )} [1 - /(log S -^_)], (89) 

(n — 1)1 

where the function f(y) is defined by Eq. ([f4]) (see also Fig. |). 

The only difference between periodic and Dirichlet boundary conditions 
is in the constant g± (|88| ) where the term < rlogr >= 2(1 — log 2) is added 
for the latter. 

As above when s — > oo the function /(logs — g±) goes to zero and the 
true asymptotics of the n-th nearest-neighbor distribution with n > 1 for all 
boundary conditions is 

„n— 1 

i.e. it coincides with the (n — l)-th nearest-neighbor distribution for the Pois- 
son distribution. A simple physical explanation of this result is the following. 

We are interested in the solutions of Eq. (0) when unperturbed levels, ej, 
are independent random variables. Among all configurations of Cj there are 
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cases where two unperturbed levels e\ and e 2 are very close to each other. In 
such a case Eq. (|7]) reduces to 2 terms 

sV^r - (91) 

which has a simple solution 

E = 6ir2 + 62ri . (92) 
n + r 2 

As we assume that the difference e\ — e 2 is very small, the value of E will 
also be very close to both unperturbed levels, e\ and e 2 . The equation fl92]) 
can be reversed and having one unperturbed level, say e\, and new level, 
E, very close to it, one can always find the position of another unperturbed 
level, e 2 , to fulfill Eq. (0) in that approximation. We shall call 2 very close 
unperturbed levels with a new level inside a dipole configuration. Now let 
us consider the probability that 2 unperturbed levels e% and e 2 are at the 
distance s with n — 1 unperturbed levels inside this interval. As unperturbed 
levels are independent this probability is given by Eq. (|90|). 

We represent in Fig. |6] a possible configuration for 2 unperturbed levels 
(short thin vertical lines) separated by a large distance s = t\ — e 2 with n — 1 
levels between them. For any such configuration there exists a configuration 
with 2 new energy levels at E 2 and E\ very close to e 2 and e\ respectively. 
This is true because, as we have pointed out, it is possible to construct two 
unperturbed energy levels e 2 and e[ (indicated by dashed lines in Fig. |(J) 
such that two pairs (e' 2 ,e 2 ), form the dipoles. For that configuration 

it is clear that E 2 and E% are two perturbed levels with n perturbed levels 
in between. When s —>■ oo it is physically clear that the other levels will not 
influence these dipole configurations which explains Eq. (p0|). 

It is clear that this reasoning cannot be applied to the nearest-neighbor 
distribution (as it requires at least 2 unperturbed levels inside the interval s) 
and the asymptotics of P(s) given by Eq. (|77|) is quite different from Eq. (|90|). 

To check the asymptotic formula (^) we compute numerically the n-th 
nearest-neighbor distributions till n = 9 for a rectangular billiard with size 
4 x 7r with periodic boundary conditions and for different coupling constants. 
First, we find the best fit of the integrated n-th nearest-neighbor distribution 
N n (s) = P n (t)dt in the form aexp(bs)s c . In Fig. [7| we plot values of b 
and c obtained by this fit. The lower curve (dots) shows the values of the 
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constant b, which as expected is the same for all values of n and is equal 
with a good precision to —1. The upper curve (squares) shows the exponent 
c which according to ( |59"D is expected to be n — 1. This is indeed the case for 
the higher values of n. The small deviations from this expected value for the 
lowest n comes from the fact that the function / has to be taken into account 
since for large s, /(logs — g±) behaves like 1/ logs. To illustrate accuracy of 
Eq. ( [39l) we present in Fig. |]the results of numerical computations for ^(s) 
and Ps(s) for a billiard with periodic boundary conditions and renormalized 
coupling constant, A' = 1. Error bars in this figure indicate statistical errors. 
It is clear that the asymptotic formula fl89|) describes well large s behavior 
ofP„(s). 

The same checks have been performed for billiards with Dirichlet bound- 



ary conditions. The results are presented in Fig. and Fig 10. Again the 



theoretical asymptotics is in a good agreement with numerical results. 

5 Conclusion 

The starting point of our investigation is Eq. (|7]) 



N 



r 



j 



, E - e j 

3=1 J 



1. (93) 



We assume that (i) all ej are independent random variables, (ii) the residues 
Tj are real positive, and we compute the n-th nearest-neighbor distribution 
P n (s) of the solutions, E, in the limit N — > oo. The exact formulas are quite 
cumbersome and we dwell on the asymptotic behavior of P n (s) at large s. 
Our main results are the following. 

The asymptotics of the nearest-neighbor distribution, P(s), is given by 
Eq. ( |77|) and has the form 

P(s) = / [<C >) =- (94) 

Here v sp is determined by the equation 

- < rEi{rv sp ) >= — (95) 
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where Ei(x) is the usual exponential integral fl63| ) and A' is the renormalized 



coupling constant (48). The notation < f(r) > indicates the mean value over 
all residues. For periodic boundary conditions < / (r) >= /(l). For Dirichlet 
conditions < f(r) > depends on the ratios of the coodinates of the scatterer 
to the corresponding sides. When these ratios are non-commensurable irra- 
tional numbers < f(r) > is defined in (|SlD. When they are rational numbers 
< f(r) > should be computed as in ([82]). 

The n-th nearest-neighbor distribution with n > 1 when s — > oo has the 
following asymptotics 

P n (s) = (^Jje- (96) 

which depends neither on the residues nor on the boundary conditions. For 
finite values of s there are slowly decreasing corrections to these formulas 
indicated in Eqs. (ff5[) and ( p9f) which are important for accurate comparison 
with results of numerical calculations. 



The above results together with the results of Ref. |L3| prove that spectral 
statistics of generic rectangular billiards with a small-size scattering center 
inside is of special (intermediate) type characterized by 2 important prop- 
erties: (i) level repulsion at small s and (ii) exponential decrease of the 
nearest-neighbor distribution at large s. For the authors' knowledge this is 
the first example of a dynamical system where the intermediate character of 
the spectral statistics can be proved rigorously. 

Appendix 

The purpose of this Appendix is the computation of two main integrals 

oo 

Ji(ai,a 2 ) = ( / + /)[l-exp(i— i — )]de, (97) 



-oo 1 



and 



J (ai,a 2 )= / [1 — exp(i- i — )]de. (98) 

Jo 1 — e e 
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As in [|TJ| we first find the difference of the derivatives over a,\ and a 2 

oo 

(dx - d 2 )Ji(a 1 ,a 2 ) = -i( / + / )- TJ r exp(z'-— ^ z— ), (99) 

V J e(l-e) 1-e e 

— oo 1 

where denotes the partial derivative with respect to the aij. Changing the 
variable e = t/(l + 1) one gets 

(di - d 2 )J l {a l , a 2 ) = e^-^K « 2 ), (100) 



where 



ii(oti, 012) = i I — exp(— ia\t + i— ). (101) 



The last integral is well defined for complex a when Im (ai) < and 
Im (a 2 ) > 0. In this domain (see |23| ) 



<M«i, a 2 ) = 2iK (2^^T 2 ) (102) 

where K (x) is the modified Bessel function of the third kind. At real a the 
limiting function is discontinuous and depends on the sign of a 

-Ti (2^/— «ia 2 ) when a\ < 0, a 2 > 

i 2iK (2^/aia 2 ) when a x > 0, a 2 > , , 

<p l (a 1 ,a 2 ) = { rr(2) fo X , ~ _ • (103) 

ttHq (2y/ —a>\a 2 ) when a,\ > 0, a 2 < 

2iK Q (2y / ai02) when ai < 0, a 2 < 
The function Jo(a) obeys the similar equation 

(ft - ft) J (ai, a 2 ) = e^- a2 Vo(ai, « 2 ), (104) 

where 

0o ("l, oi 2 ) = -i I exp(i«it - i^) = -2iK (2 v /«ia 2 ), (105) 
Jot t 

which is defined in the region Im (a-i) > and Im (ct 2 ) < 0. At real a, 
(f>o(ot) = 4>l(a) or explicitly 

-ttHq {2y/— «i« 2 ) when ai < 0, a 2 > 

= ■ -2^0 (2 y^) when a x > 0, a 2 > (10g) 

7riifQ (2-y/— atia 2 ) when ai > 0, a 2 < 

— 2ii^ (2^/aia 2 ) when < 0, a 2 < 
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We note also the expression for the sum of 0<r(a) 

0i («i, a 2 ) + 0o(«i, oc 2 ) = vr(sgn(«i) - sgn(a 2 ))^o(2v / -«ia2), (107) 

where sgn(x) denotes the sign of x. 

The knowledge of permits to write down a linear partial derivative 

equation for J a (a) 

(di - d 2 )Ja(ai, a 2 ) = Z a (ai, a 2 ) (108) 

where 

Z a (a x , a 2 ) = e^-^K, a 2 ). (109) 
The general solution of this equation has the form 

J a (ax, a 2 ) — J a (0, a x + a 2 ) + / Z a (t, a x + a 2 - t)dt. (110) 

The initial values J CT (0, a) can be computed directly from the definitions (|97|) 
and (^). By changing the variable etot = — 1/e one gets 

/0 rco poo i, 

+ J )(1 - e^de =j e lQi ), (111) 

where ^ is the principal value of the integral. Similarly 

f 1 f °° rlf 

J o (0,a) = jf (l-e- iQ / e )rfe = jf _(l- e -**). (112) 
The final expressions for J a (a) are the following 

/oo j± pai 
_ (1 _ e ,(a 1+Q2 )t )+ / gi^-ai-^^^^ + ^-t)^ 

(113) 

and 

/oo j. /-ai 
_(1 _ e -(«!+« 2 )*) + y e ,(2t- Ql -a 2 ) 0o(t; ai + ^ _ t)dt 

(114) 
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According to Eqs. (j9S|) and (TfjTfy the functions Jo(a) and Ji(a) considered as 
functions of complex ot\ and «2 are analytical functions in different regions: 
Im («i) > and Im (a 2 ) < for the former and Im (o^) < and Im (c^) > 
for the latter. They can be continued to complex planes with the cuts 
represented in Fig. |l|a and |I|b for J\ and in Fig. |l]c and |l|d for J . At real 
values of a\ and a 2 they have a discontinuity along the axis a.\ — and 
a 2 = 0. 

The sum Jo(a) + Ji(a) obeys an equation similar to (|100|) where instead 
of function ^i(a) one substitutes the sum (|107|) 

(d 1 -d 2 )[J 1 (a 1 ,a 2 ) + J («i,« 2 )] = ^-^Z{a u a 2 ) (115) 

where 



Z(a u a 2 ) = 7r(sgn(«i) - sgn(a 2 )) J (2 v / -«i«2)- (116) 

The solution of these equations can be done exactly as above and it can be 
represented as a sum of two discontinuous functions 

J (ai, a 2 ) + Ji(«i, a 2 ) = sgn(cti)i?(ai, a 2 ) + sgn(a 2 )-R*(a2, «i), (117) 
where 

R(a h a 2 ) = n [ 1 J (2y/-t(ai + a 2 - t))e l{2t - ai - a2) dt. (118) 
Jo 

One can check that this expression coincides with Eq. (A24) of ||13|| . 

The following useful symmetry properties can be checked directly from 
the definitions (0)-(0) 

J a {a 2 , ax) = J a (-a u -a 2 ) = J*(ai, a 2 ), 

4> a {a 2 , ai) = CT (-ai, -at 2 ) = "2), (119) 

R(—o>i,—a 2 ) = —R*(o>i,a 2 ). 

For further references we present the behavior of J a (a) and <f) a (oc) at small 
real a 

Ji(a 1} a 2 ) = |(|a 1 | + \a 2 \) (120) 
+w*i(l - 7 - log |oti|) - ia 2 (l - 7 - log \a 2 \), 

7T 

J (oci,a 2 ) = -(|«i| + \a 2 \) (121) 
-ia x (l - 7 - log |ai|) + ia 2 (l - 7 - log \a 2 \), 
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and 

7T 



>l 



Oil, ol 2 ) = -(sgn(ai) - sgn(a 2 )) - 2(27 + log |«i| + log \a 2 \), (122) 



7T 



<f>o(ati,a 2 ) = -(sgn(ai) - sgn(a 2 )) + i(2?f + log |ai| +log|a 2 |), (123) 

where 7 is the Euler constant. 

Due to the above mentioned analytical properties these expressions 
(though they have apparent discontinuities) can be rewritten as following 
analytical functions 

Ji(ai, a 2 ) = zai(l - 7 - log(zai)) - ia 2 {\ - 7 - log(-za 2 )), (124) 

J (ai, oc 2 ) = -z«i(l - 7 - log(— ieti)) + ia 2 (l - 7 - log(za 2 )), (125) 

and 

0i (ai, a 2 ) = -1(27 + log(iai) + log(-ia 2 )), (126) 

0o(«i, a 2 ) = ^(27 + log(-iai) + log(ia 2 )). (127) 

where we assume the usual definition of the logarithmic function with a cut 
along the real negative axis. 

Similarly as it was done in [I3"| for the sum Ji(a) + Jo (a) one can also 



obtain higher order terms of the expansions of J a (a) in power of a. Both 
functions J\ and J given by Eqs. (|113|) , ( |114j) contain terms proportional to 
log(«i + a 2 ) but one can check by direct series expansions that in the corre- 
sponding sums these terms all cancel and the only logarithmic contributions 
to these functions are the ones presented in Eqs. ( |120| ) and ( |121| ) as it should 
be to ensure the analytical properties of J a (a). 
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Figure 6: Schematic representation of the dominant contribution to the n- 
th nearest-neighbor distribution. Thin vertical lines are unperturbed energy 
levels. Thick lines represent new energy levels. 2 vertical dashed lines indicate 
2 unperturbed energy levels which are added to construct encircled dipole 
configurations. 
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Figure 7: Values of b (dots) and c (squares) in a fit of N n (s) for periodic 
boundary conditions under the form aexp(bs)s c . 
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Figure 8: i^s) (the left curve) and Ps(s) (the right curve) for a billiard with 
periodic boundary conditions. 
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Figure 10: p2(s) (the left curve) and Ps(s) (the right curve) for a billiard 
with Dirichlet boundary condition 
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